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ABSTRACT 

This paper compares reported dynamic analyses for evaluating the 
steady-state response and stability of free-piston Stirling engine 
(FPSE) systems. Various analytical approaches are discussed to 
provide guidance on their salient features. Recommendations are 
made in the recommendations remarks for an approach which 
captures most of the inherent properties of the engine. Such an 
approach has the potential for yielding results which will closely 
match practical FPSE-load systems. 

I. INTRODUCTION 

Historically, the dynamic analyses of free-piston Stirling engine 
(FPSE) systems have focussed on two areas. Most of the analyses 
have emphasized the use of Laplace transformation and classical 
control techniques to determine the engine operating frequency and 
other dynamic parameters, notably, the piston-displacer oscillation 
amplitudes and their relative phase angles [1,2,3]. The other area of 
activity has been the determination of the conditions for engine 
stability by means of Laplace transformation [4] or the state-space 
technique [5,6]. In one instance [7], the operating frequency, piston 
and displacer amplitude, and other thermodynamic state variables are 
found by a linear harmonic analysis (LHA). This technique 
represents periodic variables with harmonic functions. The underlying 
reason for the bulk of the dynamic analyses is the prediction of the 
FPSE performance. The foregoing activities are essential to the 
design of a particular Stirling engine, and can enhance one s 
understanding of the behavior of FPSE’s. 

'Diis paper reviews the existing literature on the dynamic analysis of 
FPSE systems. The purpose is to discuss the various analytical 
methods used, and provide guidance on their relative merits. The 
recommendations state the approach and salient features of the 
FPSE-load systems which will yield the most practical results, in a 
most expeditious manner. 

II. ANALYTICAL FORMULATIONS 

The discussions here pertain to the underlying assumptions and the 
equations of motion used for the dynamic analysis of FPSE s. For the 
papers reviewed and in this publication, the word system denotes a 
FPSE connected to a load. Where possible, the nature of the load is 
identified. 

II-l. FPSE SYSTEMS 

AH the analyses found in the available literature on ITSE dynamics 


emphasize the engine, with only allusions to its connected load. In 
some cases, the system consists of a single cylinder engine with a 
connected but unspecified load [1 ,5,6,7] modeled as a dashpot. In 
other instances, the system is a FPSE driving a linear alternator [2,3]. 
In reference [2], a connected load is implied but not explicitly stated. 
If the sole desire is the calculation of the engine dynamic parameters 
which impact the thermodynamic analysis, the system is usually 
confined to the engine itself [4]. 

Generally, the dynamic analysis of the FPSE system is expedited by 
assumptions which facilitate the solution of the system equations. 

II-2. ASSUMPTIONS 

The following simplifying assumptions are commonly found in the 
analyses: 

(1) . Schmidt’s thermodynamic analysis is assumed. Hence, (a) the 
piston and displacer motions are sinusoidal, as is the resulting 
working space pressure; (b) the working fluid obeys the ideal gas law 
and undergoes isothermal expansion and compression. 

(2) . The system can be studied via linear analysis methods since (a) 
the working gas behaves like a linear spring and, (b) the connected 
load is assumed to be nearly linear. 

(3) . The gas pressure in the working space (expansion and 

compression spaces, heater, regenerator and cooler volumes - Fig 1) 
is spatially constant but time variant. 

(4) . The average working space pressure is balanced by the average 
pressure of the bounce space and other gas springs. Thus, the 
average force on the power piston, displacer and cylinder casing are 
zero. Consequently, the average positions of these elements are 
stationary. 

II-3. EQUATIONS OF MOTION 

A schematic of a single cylinder FPSE used by Rauch [1] is shown in 
Fig. 1. The dynamic equivalent model is illustrated in Fig. 2 . The 
FPSE is represented as two masses, namely, the displacer (subscript 
’D’) and the power piston (subscript ’P’), which are coupled by 
gaseous springs (denoted ’K.’) and dampers (denoted C). Unlike Fig. 
2 in reference [1], an additional damping coefficient which couples the 
piston and displacer in the compression space is included here for 
completeness. Positive displacements of the piston and displacer are 
considered upward motions from their static equilibrium positions in 
Fig. 2. 


The dynamics must be self-excited via the thermodynamics so as to 
induce the engine operation. That is, the dynamics of the piston and 
displacer motions must generate a pressure force which maintains 
their steady, periodic motions. 


[K ? ] 
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The piston, displacer and casing arc three key elements of the FPSE. 
A sinusoidal steady oscillation is commonly assumed for these 
elements. Urieli and Berchowitz [9] show that operating the FPSE at 
a frequency in excess of the natural frequency of any of its elements 
will result in a reduced amplitude of oscillation for that element. One 
mode of FPSE operation is the removal of power out of the casing. 
For this mode the operating frequency is much greater than the 
natural frequency of the piston. Another mode is the removal of 
power from the piston. This mode is common for electrical output 
applications in which an alternator "load" is installed between the 
piston and casing. In this case, the operating frequency is greater 
than the natural frequency of the casing. This mode of operation is 
preferred for space power applications (and many others), due to the 
simplification obtained by a relatively heavier and, hence, stationary 
casing. 

Based on the foregoing discussions, only two degrees of freedom are 
required to describe the dynamics of the FPSE. Hence, regarding Fig. 
2, the general equations of motion are summarized in Eq. (1): 

[M][XMC][XMK][X]«IF(t)] (1) 


where and K pp p are the total thermodynamic stiffness 

coefficients of the displacer and piston, respectively. The term Kppj 
is the thermodynamic stiffness on the displacer due to the piston 
motion. Similarly, Kp£ T is the thermodynamic stiffness on the piston 
due to the displacer motion. 

The terms in [Kj] are functions of engine geometry, operating mean 
pressure and the expansion to compression space temperature ratio, 
T e /T c . Generally, these stiffness coefficients are uniquely different 
from each other. In particular, the coupling terms K DPT and K PDT 
are not equal, making [Kj] asymmetric. This asymmetry is the source 
of self-excitation in the practical engine. However, when the working 
space temperatures are equal, that is T e /T c is unity, the coupling 
terms K DPT and K PP j are equal and [K T ] becomes symmetric. 

The thermodynamic stiffness matrix [K T ] may be split into an 
isothermal term [KjjJ and a temperature ratio dependent term 
[K Tt ]. Thus, Eq. (3) becomes: 

[F T (t)]= -[Cj][X]-[K Ti ♦ K xt ][X] (5) 


l>ie matrices [M],[C] and [K] represent the system masses, damping 
and stiffness coefficients external to the thermodynamic cycle. The 
forcing vector [F(t)J is the sum of the forces due to thermodynamic 
cycle pressure, [F T (t)] f and appropriate time dependent external 
forces, [Fp(t)), on the piston and displacer. This is shown in Eq 2: 

[FcohFtUMFeCO] (2) 


The thermodynamic forces can be represented as functions of the 
state variables [X] and [X], which represent the displacements and 
their time derivatives, respectively, for the piston and displacer, as 
shown in Eq. 3: 

[F T (t)]-|C T ]|Xj-[K T ](X] (3) 

where [C T ] and [K T ] are, respectively, the thermodynamic damping 
and stiffness matrices which may be nonlinear. Negative signs 
associated with [K T ] and [C T ] indicate the restoring nature of [F T (t)J. 

l>ie matrix [C-p] is, in most analyses, given only passing attention. It 
represents the heat exchangers flow losses which, in many cases, are 
the dominant engine losses. Further, it may be shown that this 
damping is a non-linearity which stabilizes the amplitude of oscillation 

HI- 

The [K t ] matrix represents the dominant pressure forces acting on 
the piston and displacer. It may be written in terms of the piston 
area (Ap), displacer rod area (A R ) and the partials of pressure with 
respect to positions, as: 


The [KjJ matrix is symmetric and does not contribute to the self- 
excitation. The [K Tt ] matrix is asymmetric, retains the temperature 
ratio dependence, and is the primary contributor to the self-excitation. 

There are several ways of treating the thermodynamic force [F T (t)]. 
One approach is to combine the thermodynamic matrices [K T ] and 
[C-j] with their external counterparts [K] and [C] on the left hand side 
(LHS) of Eq. (1), which becomes: 

lMJ[X]*{C*Cj][X]+[K+K T JpC]«[F£(t)) (6) 

A second approach is to model the thermodynamic cycle as a forcing 
term in the right hand side (RHS) of Eq. (1). Combining Eqs. (1) 
and (2) yields Eq. (7): 

(M][XMC][XHK][XHF T (t)MF E (t)] (7) 

In this case Eq. (3) must also be satisfied for the complete solution. 

A third method is to split the thermodynamic model, in which case 
[C T ] and [K Ti ] terms of Eq. (5) are combined with [C] and [K], 
respectively, on the LHS of Eq. (1). The [K Tt ] term remains on the 
RHS of (1) as a forcing term. Eq. (1) then becomes, 

[M][Xj^C-C T ][X] + [K*K Ti HX]-[F pI (t)]*[F E (t)] (8) 

where the force representing the relation between the thermodynamic 
cycle and piston and displacer motions, Eq. (9), must also be satisfied. 
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IF p t( 0] - -[K Tt ][X] (9) 


(K^K t ] = 


(13) 


Various forms of Eqs. (6), (7) or (8) are used in the literature. The 
equation used and unique assumptions of each reference are 
discussed below. 


Rauch [1] splits the thermodynamic model and, thus, uses Eq. (8). 
The cylinder is assumed to be stationary. The engine modeled has no 
external springs. This implies that the elements in the expanded 
stiffness matrix of Eq. (10) are solely the isothermal counterparts of 

[Kill- 


[K+K Ti ] = 


“Kc * K P +K C 


( 10 ) 


The common, off-diagonal term represents the displacer-piston 
coupling spring stiffness. The terms Kq and Kp, respectively, denote 
the displacer and piston springs to ground. Note that the diagonal 
terms are the sum of K D and K^ for the displacer, and Kp and Kq for 
the piston. 


Reference [1] correctly points out that the elements of the total 
damping matrix [C + C T ] are dependent on both the engine-connected 
load and internal windage. The damping dissipates the energy input 
to the piston and displacer. Maximization of system efficiency 
requires minimization of damping other than that due to the load. 
Reference [1] embeds the effects of connected load in the piston 
damping coefficient, C p . Thus the damping matrix may be expanded 
to Eq. (11). 


[OC T J = 


Cd+Cc.-Q; 

-C c .Cp+Cc 


( 11 ) 


To simplify the analysis, reference [1] neglects the coupling term C c 
in the damping matrix. The forcing terms of (8) are assumed to be 
entirely due to the thermodynamics. Therefore, [F E (t)] is zero. The 
forcing vector [F pX (t)] is expressed in terms of a time dependent 
pseudo-pressure acting on the piston and displacer rod areas as shown 
in Eq. (12): 


[F pX (t)] = 



( 12 ) 


The constraint of Eq. (9) on the equations of motion requires 
equating the RHS of Eqs. (9) and (12). This is not explicitly evident 
in reference [1J. 

Unlike reference [1], Redlich and Berchowitz [2] incorporate the 
Schmidt thermodynamics into the total matrix. Thus, their equations 
of motion are similar to Eq. (6). The total stiffness matrix is 
expressed in Eq. (13), 


K d ,a p 

a T * Kp 


where K d represents both the external and thermodynamic (K DDT ) 
effects. The Kp includes only the thermodynamic (K PPT ) effects. 
The term a p is due to thermodynamic effects and represents the 
gaseous force exerted on the displacer rod area due to the piston 
motion. It is equivalent to K DP p in Eq. (4). The term a x a 
thermodynamic coupling between the displacer motion and the piston 
force, denotes K PDT in Eq. (4). Reference [2] shows that K P , K d and 
the stiffness coefficients a P and ap are functions of the displacer rod 
area, engine cylinder area, the gas pressure and the expansion and 
compression space temperatures. This observation is consistent with 
Eq. (4). Thus, values of these parameters may vary for Stirling 
engines with different geometric configurations and operating 
temperatures and pressures. 


The damping matrix for reference [2] is shown in Eq. (14): 


[OC T ] = 


D d . 

0, 



(14) 


The damping coefficient, D d , of Eq. (14) embodies all viscous forces 
on the moving gas. The term, D p , includes the effect of any piston- 
connected load. Both D d and include "incidental irreversibilities". 
An example of this is gas spring hysteresis. There is no explicit 
damping coupling term between piston and displacer in Eq. (14). 
Reference [2] equates [F E (t)] to zero since there are no external 
forces. 

The mechanical analog representation of the Stirling engine in Fig. 2 
is nearly identical to that used by Das and Bahrami [3]. Also, the 
underlying assumptions in references [1] and [3] are similar. Hence, 
Eq. (8) represents the system equations for reference [3]. Das and 
Bahrami quote reference [1] assertions that the piston-displacer 
coupling is weak and can be neglected in a well designed and efficient 
engine. They support this by noting the relatively heavier mass of the 
load compared to that of the piston, and the relatively dominant 
stiffness of the engine gas spring action. However, Das and Bahrami 
do not explicitly ignore the coupling terms in formulating the 
equations of motion. 

The formulation of the thermodynamic stiffness matrix, [Kp], and the 
dynamic equations by Benvenuto, et al., [4] is similar to that by 
Redlich and Berchowitz [2], namely, Eq. (6) with a few exceptions. 
Reference [4] more rigorously develops the thermodynamic damping 
matrix [Cp]. This matrix is shown to be a non-linear function of the 
piston and displacer velocities. Both the direct and coupling terms 
are included in the formulation. 

Cichy, Carlini and Kucharski [5,6] use a formulation similar to Eq. 
(7). In their analysis, the [K] and [C] matrices are due entirely to 
components external to the thermodynamics. The displacer stiffness 
to ground is not explicitly shown by Carlini, et al. This is consistent 
with their dynamic representation of the Beale model 10B Stirling 
engine. Their system has no external forces. Thus, [Fp(t)] is zero. 
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Their thermodynamic forcing term, [F T (t)|, is similar to Eq. (15), 
which has been rewritten to be consistent with Fig. 2., 


[F T (t)J- 


-A R +Ap, A p P e (t) 

°, -A p J|p c (t) 


(15) 


where P e (t) and P c (t) are, respectively, the expansion and 
compression space pressures. Thus, the thermodynamic forces are 
modeled as non-homogencous terms in the dynamic equations. The 
natural frequencies calculated with this model will differ from reality, 
since the "stiffness" effects of the working space are incorporated in 
the forcing terms [F T (t)]. The constraint of Eq. (5) on the equations 
of motion is not explicitly recognized. This requires equating the 
RHS of Eqs. (5) and (15). 

The analysis by Chen and Griffin [7] is an extension of their previous 
linear harmonic analysis work [8] which models the thermodynamics 
of the Stirling cycle. In reference [7] they have coupled the mechanical 
dynamic equations of motion with the differential equations 
describing the thermodynamic processes. Thus, Eq. (6) is most 
similar to their analysts. The major difference in their formulation is 
that the order of the matrices [X], [K T ] and [C T ] in Eq. (3) are 
expanded to explicitly include thermodynamic and dynamic variables, 
as well as losses and adiabatic effects in the thermodynamic cycle. 

Literature search reveals that the selected approach for solving the 
equations of motion depends on the objective of the author. The 
commonly used methods are discussed next. 


Dynamic stability of the FPSE is obtained from the system 
characteristic equation. Application of Laplace transformation and 
Routh’s stability criterion gives the impact of parameter variations on 
engine stability. 

Benvenuto, et al. [4] approach is to linearize the working gas pressure, 
pressures in the gas springs and pressure drop in the heat exchangers. 
The result is a system of linear and homogeneous equations of 
motion with constant coefficients. These equations are Laplace 
transformed into the complex domain. Solution of the polynomial 
characteristic equation yields the conditions for stable dynamic 
behavior of the FPSE. 

Cichy and Carlini [5] and Carlini and Kucharski [6] employ the state- 
space formulation of the equations of motion for their analyses. 
Reference [5] uses a two-part approach for determining the 
thermodynamic parameters. Initially, a spatially uniform sinusoidal 
pressure variation is assumed for the working fluid. This is the so- 
called "zero-approximation" analysis. The pressure losses between the 
compression and expansion volume are included in a subsequent 
analysis. This is termed the "first approximation" analysis. For either 
type of analysis, the thermodynamic parameters and transient 
response are obtained for a unity step input and sinusoidal pressure 
input. 

Carlini and Kucharski [6] apply the modal transformation and 
eigenvalue analysis to obtain the control requirements for, and 
parametric effects on the engine dynamic behavior, following both a 
step function and harmonic (or sinusoidal) input types of excitation. 


III. METHODS OF SOLUTION 

Analysts have their preferences for the technique employed in solving 
the dynamic equations. Time-domain analysis, either by direct 
solution of the differential equations or matrix formulation of the 
equations, is commonly used to evaluate the operating frequency, the 
displacer-piston displacements and relative phase angles. Other 
analysts opt for a combination of time- and frequency-domain 
analyses in establishing steady-state system stability and/or criterion 
for its occurrence. 


Chen and Griffin [7] assume that all system variables are harmonic in 
nature. They then reduce the differential equations to a system of 18 
homogeneous algebraic equations. These are solved for the operating 
frequency by means of a secant-bisection root-finding method. This 
method ensures stable convergence in a minimum time. The 
homogeneous set of equations is converted to non-homogeneous type 
by assuming a value for power piston amplitude. The equation 
representing the piston displacement is deleted in the set of 
homogeneous equations. 


Rauch [1] obtains the amplitudes of the piston and the displacerand The influence of the piston motion on the remaining equations is 

their relative phase angles by evoking the Schmidt sinusoidal motion replaced by a forcing term which is proportional to the piston dis- 

for the exerted force, and the piston and displacer responses. placement. Standard matrix algebra is used to solve for 17 of the 18 

Substitution of the force and responses into the matrix equations of system variables, based on an assumed piston amplitude. A unique 

motion, and subsequent application of Gaussian elimination yields the solution for the piston amplitude is found by matching the power 

solution for the basic equations at the desired frequency. produced by the engine to the load power requirement. The final 

solutions are obtained by successive iteration of the above solution 
Rcdlich and Berchowitz [2] use Taylor series expansion to linearize techniques, 
the equations of motion about a steady-state operating point. The 

resulting equations are Laplace transformed. Nyquist criterion is IV. CONCLUDING REMARKS 

applied to the characteristic equation to evaluate the necessary 

condition for system stability, and criteria for engine start-up, and The primary objective of the reviewed dynamic analyses is to 
maintenance of piston and displacer resonances. determine the dynamic behavior of FPSE systems. Several authors 

formulate their analyses such that system stability and parametric 
Das and Bahrami [3] assume a sinusoidal forcing function and, hence, constraints for stable operation may be evaluated. However, with few 

displacer and piston responses. These arc substituted into the exceptions, stability has been defined in the classic sense, i.c., with 

equations of motion which are rewritten in complex exponential form respect to static equilibrium, rather than "periodic" stability in which 

to determine the piston and displacer steady-state responses. the system converges to a stable limit cycle. 
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Frequency domain analysis in the form of Laplace Transformation, 
with inverse transformation for time-domain response, has been 
effective in computing the basic thermodynamic parameters. 

V. RECOMMENDATIONS 

The dynamic behavior and stability of FPSEs depend on the total 
system, including the thermodynamics and load. Therefore, a 
complete analysis should incorporate detailed representations of both 
the thermodynamics, dynamics and load. In particular, for engines 
with a connected linear alternator, the model should include the 
electrical dynamics of the alternator and electrical load. Also, the 
thermodynamic model should account for both adiabatic effects and 
major losses directly affecting the working space pressure amplitude 
and phase. 

The thermodynamic stiffness and damping coefficients, including the 
coupling terms, should be included in the homogeneous equation to 
yield a more accurate prediction of the system behavior. Thus, Eq (6) 
is recommended since all first-order dynamic effects are satisfied. 
Further, the use of Eqs (3) and (4), with non-linearities included in 
the matrices, for describing the thermodynamic system, is deemed 
sufficient. However, including the detailed thermodynamic equations 
with the dynamic analysis [7] is appropriate where thermodynamic 
performance is the primary concern. 

The major source of self-excitation is contained in the stiffness matrix, 
while the major stabilizing influence is contained in the damping 
matrix. Both matrices can be shown to contain non-linear 
coefficients. Better understanding of the nature of various non-linear 
effects and their influence on system behavior is needed. Inclusion 
of the non-linearities should produce more realistic results. Although 
approximate, linearized solutions to the non-linear equations of 
motion can yield useful results, detailed solutions will require 
numerical analysis. 

Frequency-domain analysis, using a state-space technique, facilitates 
a detailed sensitivity analysis of the effects of parameter variations on 
the engine dynamic behavior. This is useful, particularly during the 
engine design stage, in establishing acceptable stability margins. 

Only a few of the analyses [2,7,9] have been compared to 
experimental results. Experimental validation of the dynamic analysis 
is recommended to enhance the analytic models and increase 
confidence in their predictive capability. 

In contrast to dynamic stability, the transient response of the FPSE 
system will need to be evaluated by use of time domain analysis. 
Transient response, a potential operating mode of FPSE, has not yet 
received wide-spread attention. 
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GLOSSARY: 

Symbols: 

A - Area 

C, D - Damping coefficient 

FPSE - Free-Piston Stirling Engine 

K - Stiffness coefficient 

P - Pressure 

Matrices /Vectors: 

[C] - Damping matrix 

[F(t)] - Forcing vector 

[K] - Stiffness matrix 

[M] - Mass matrix 

[X],[X],[X] - State variable vectors: position, velocity', acceleration 
Subscripts: 

C - Displacer - piston coupling 

c - Compression space 

D, d - Displacer 

DDT - Total thermodynamic, for displacer stiffness 

DPT . = dp = Thermodynamic stiffness on displacer due to 

piston motion 
E - External 

e - Expansion space 

p - Piston 

PDT - s ay = Thermodynamic stiffness on piston due to 

displacer motion 

PPT - Total thermodynamic, for piston stiffness 
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